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Abstract 


The  Lax-Ricthmyer  theory  is  used  to  study  the  error  amplification 
properties  of  18  space  marching  finite  difference  schemes,  for  the  1-D 
nonlinear  inverse  heat  conduction  problem.  A non-dimensional  param- 
eter Q.  involving  the  time  step  AT  the  effective  thermal  diffusivity  a, 
and  the  distance  I from  the  sensor  to  the  active  surface,  is  found  to 
provide  a measure  of  the  numerical  difficulty  of  the  inverse  calcula- 
tion. All  18  schemes  are  unstable  and  blow-up  like  10^^,  where  the 
constant  A depends  on  the  particular  numerical  method.  There  are 
substantial  differences  in  the  A's  however,  and  some  newly  constructed 
algorithms,  employing  forward  time  differences  at  non-adjacent  mesh 
points,  are  shown  to  produce  relatively  low  values  of  A.  Using  synthetic 
noisy  data,  a nonlinear  reconstruction  problem  is  considered  for  which 
Q = 25.  This  problem  simulates  heat  transfer  in  gun  barrels  when  a 
shell  is  fired.  It  is  shown  that  while  most  of  the  18  schemes  cannot 
recover  the  thermal  pulses  at  the  gun  tube  wall,  two  of  the  new  meth- 
ods provide  reasonably  accurate  results.  A tendency  to  underestimate 
peak  values  in  fast,  narrow  thermal  pulses  is  also  noted. 

Key  Words:  nonlinear  heat  flow,  marching  difference  schemes,  er- 
ror amplification,  numerical  e.xperiments. 

AMS  (MOS)  subject  classifications:  35R25,  65M30. 
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1 Introduction 


The  inverse  heat  conduction  problem,  (IHCP),  whereby  surface  temperature 
and  gradient  histories  are  calculated  from  histories  measured  at  an  accessible 
interior  location,  remains  a basic  problem  in  many  areas  of  heat  transfer. 
On  page  2 of  the  monograph  by  Beck.  Blackwell,  and  St.  Clair,  [1],  the 
authors  indicate  a variety  of  apphcations,  and  they  estimate  that  some  300 
papers  have  been  written  in  this  general  area  as  of  1985.  A substantial 
mathematical  literature  has  also  developed  around  the  IHCP  during  the  last 
30  years  or  so.  as  may  be  seen  from  the  bibhography  in  Hao  and  Gorenflo. 
[6].  An  interesting  new  approach  based  on  systems  theory,  together  with 
a brief  overview  of  current  methods  and  some  further  references,  may  be 
found  in  hlarquardt  and  Auracher.  [9]. 

The  present  paper  focuses  exclusively  on  space  marching  finite  difference 
methods  for  the  one  dimensional  problem  in  a medium  with  temperature- 
dependent  thermal  properties.  Several  such  methods  are  presented  in  Beck 
et  al,  [1].  Here,  we  do  not  discuss  procedures  based  on  replacing  the  heat 
conduction  ec[uation  with  an  approximating  hyperbohc  equation,  [5],  [17], 
nor  methods  combining  space  marching  with  time  marching,  [11].  Using 
the  Lax-Ricthmyer  theory,  [12],  [13],  we  examine  the  norms  of  discrete 
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solution  operators  associated  with  several  marching  schemes,  in  the  Fourier 
transform  domain.  This  enables  us  to  demonstrate  the  considerable  dif- 
ferences which  exist  between  various  methods  with  regard  to  error  ampli- 
fication, and  leads  us  to  construct  some  new  schemes  that  are  relatively 
well-behaved.  Numerical  reconstruction  experiments  with  synthetic  noisy 
data  are  then  used  to  illustrate  the  capabilities  of  the  new  schemes.  By 
comparing  with  known  ‘exact'  solutions,  we  demonstrate  credible  approxi- 
mate solutions  at  parameter  values  that  are  not  tractable  by  more  standard 
marching  algorithms. 

It  develops  that  the  role  of  time  differencing  is  paramount,  typically  over- 
shadowing the  influence  of  space  differencing,  that  some  0(At)  methods  are 
preferable  to  some  0{At^)  methods,  and  that  a certain  kind  of  0(Af)  for- 
ward time  differencing  can  dramatically  reduce  error  amphfication.  The 
beneficial  aspects  of  using  future  temperatures  have  been  noted  many  times 
in  the  literature.  However,  w'e  provide  some  quantitative  facts  that  are  inde- 
pendent of  the  particular  surface  profiles  to  be  reconstructed,  but  pertain  to 
the  schemes  themselves.  Computational  exploration  of  the  characteristics  of 
18  different  methods,  albeit  in  hmited  ranges  of  parameter  values,  indicates 
the  following.  In  these  methods,  the  mcLximum  amplification  factor  behaves 
roughly  according  to  Amax  ~ 10^^,  where  fl  = /(aAt)“C2^  / being  the  dis- 
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tance  between  the  internal  sensor  and  the  active  surface,  and  a the  effective 
thermal  diffusivity.  The  non  dimensional  quantity  Q reflects  the  underly- 
ing physical  situation,  including  the  desired  resolution  At.  and  provides  a 
measure  of  the  difficulty  of  the  particular  IHCP:  Q.  has  the  same  value  in 
each  marching  scheme.  On  the  other  hand.  A depends  primarily  on  the  time 
differencing  option  associated  with  a particular  method  and  is  essentially 
independent  of  the  physical  parameters  and  mesh  sizes.  We  show  that  some 
forms  of  forward  differencing  are  more  effective  than  others  in  producing 
lower  values  for  A.  while  one  particular  future  temperature  method  leads  to 
a larger  A and  makes  things  worse. 

Marching  difference  schemes  have  been  a mainstay  in  the  numerical  com- 
putation of  well-posed  imtidil  value  problems.  In  many  cases,  difficult  nonhn- 
ear  problems  can  be  solved  effectively  by  taking  sufficiently  small  marching 
steps  while  lagging  the  nonlinearity  at  the  previous  step.  A Courant  con- 
dition. hnking  the  sizes  of  the  space  and  time  increments,  may  need  to 
be  obeyed  to  maintain  computational  stabihty.  However,  for  ill-posed  ini- 
tial value  problems  such  as  the  IHCP.  all  consistent  marching  difference 
schemes  are  necessarily  unconditionally  unstable.  In  that  case,  no  Courant 
condition  exists  that  can  prevent  drastic  error  amplification  if  the  mesh  is 
made  sufficiently  fine,  [12,  p.  59].  The  salient  observation  of  the  present 
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paper  is  that  such  marching  algorithms  can  ‘blow-up’  at  widely  different 
rates  as  At,  Ax  J.  0.  Given  a fixed  fine  mesh,  error  amplihcation  may  be 
totally  overwhelming  in  some  methods,  while  being  sufficiently  mild  in  other 
methods  to  permit  reasonably  accurate  answers. 

Our  frame  of  reference  for  these  comparisons  is  an  artificial  but  repre- 
sentative problem  with  units.  The  object  is  to  estimate  temperature  and 
flux  histories  on  the  surface  of  an  SAE  4340  steel  plate,  using  measurements 
from  an  internal  sensor  located  0.50  mm  away  from  the  active  surface,  in 
the  presence  of  0.1%  added  noise.  The  effective  thermal  diffusivity  is  on  the 
order  of  0.004  mm^  miUisec~^ , while  the  surface  histories  are  pulse  wave- 
forms with  a rise  time  of  10  millisec  or  less.  To  accommodate  possibly 
rapid  changes,  1000  equispaced  mesh  points  are  placed  on  the  time  interval 
0 < t < 100  millisec.  ( 11  ~ 25  in  the  present  case).  This  choice  of  process 
parameters  and  time  scales  is  motivated  by  an  interesting  physical  problem 
involving  nonlinear  inverse  heat  transfer  in  gun  barrels,  and  the  pulse  wav’e- 
forms  at  the  plate  surface  are  chosen  to  simulate  thermal  histories  at  the 
inside  wall  of  a cannon  when  a shell  is  fired,  [4],  [16].  Reconstruction  of  such 
waU  histories  is  important  in  the  study  of  gun  tube  erosion,  [2].  Here,  fabri- 
cated surface  temperature  histories,  together  with  thermal  diffusivity  data 
for  SAE  4340  steel,  are  used  to  create  synthetic  histories  at  the  sensor  loca- 


6 


tion  by  numerically  solving  a direct  nonlinear  heat  flow  problem.  Uniformly 
distributed  random  noise  is  then  added  to  these  data,  and  the  performance 
of  various  marching  procedures  is  evaluated  on  this  data  set. 

2 The  direct  nonlinear  problem 

The  synthetic  data  are  generated  as  follows.  With  ~(n)  = p{u)  Cp{u) 
the  density-speciflc  heat  product.  k{u)  the  thermal  conductivity.  a{u)  = 
k{u)/^'{u)  the  thermal  diffusivity.  and  f{t)  = 300(  1 — 0.02f  - 0.12  x 10  ). 

consider  the  following  direct  problem  for  the  heat  conduction  equation: 

^■(u)  lit  = {k{ii)  0 < .r  < 1mm.  0 < t < lOOmillisec. 

u{x,Q)  = 300° K,  0 < X < 1mm. 

uil.t)  = {f{t)  + 3o00e~^'^°ksin{7ct/40)}°K.  0 <t  < lOOmiUisec. 

u{0,t)  = {300  4- 250.sin(7rt/100)}°A'.  0 < ^ < lOOmU/fsec.  (1) 

The  boundary  x = 1 is  the  active  surface,  and  u{l.t)  is  pulse  shaped,  in- 
creasing from  300°/i*  to  near  lOOO^A'  during  the  first  ten  or  so  milliseconds. 
See  Fig.  1.  A boundary  condition  at  x = 0 is  imposed  in  order  to  close  the 
system.  (The  choice  for  nfO.f)  is  artificial  and  does  not  reflect  any  particu- 
lar physical  process).  We  assume  a constant  w’  in  the  present  hypothetical 
experiment,  although  nonlinearity  in  u;  poses  no  additional  difficulties  and 


is  provided  for  in  the  marching  schemes  of  Section  3.  On  dividing  both 
sides  of  the  heat  equation  by  cj,  k{u)  is  replaced  by  a{iL)  on  the  right  hand 
side.  Using  thermophysical  data  for  SAE  4340  steel,  [14],  [15,  p.  395],  we 
postulate  the  following  linear  approximation  for  a{u)  in  gun  tubes: 

a(  )=  {0.01  — 8.0  X 10~*^(  — 250)}  mm^  milli.sec~^ , (2) 

in  the  range  300° A'  < u < 1000° A’.  In  that  range.  a{u]  undergoes  a 240% 
change.  Numerical  computation  of  this  direct  problem  was  accomplished  us- 
ing PDECOL,  [8].  This  software  package  uses  piecewise  polynomial  colloca- 
tion methods  for  spatial  discretization,  together  with  adaptive  step  selection 
for  marching  the  solution  forward  in  time.  Note  that  since  both  boundary 
temperatures  are  prescribed  in  (1),  the  boundary  fluxes  must  be  determined 
by  solving  the  direct  problem.  In  fact,  temperature  and  flux  histories  at  1000 
equispaced  mesh  points  on  the  time  interval  0 < t < Tmax  — lOOmillisec, 
were  obtained  at  several  interior  locations  Xj,  as  well  as  at  the  two  bound- 
aries X = 0 and  x = 1.  (Data  at  x = 0,  0.5,  and  1,  are  shown  in  Fig.  1). 
Each  interior  data  value  h{tk)  was  then  perturbed,  by  adding  to  it  a ran- 
dom number  drawn  from  a uniform  distribution  in  the  range  ±0.001h(i a.-)- 
In  the  inverse  problem  calculation,  interior  data  at  any  one  such  x,  may  be 
used  as  initial  values  in  space  marching  schemes  for  reconstructing  surface 
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histories.  (Some  schemes  use  u(x,.  t)  and  u{Xt  - A.T.f).  rather  than  u(x,,  #) 
and  iZx(x,.  t)).  In  a practical  setting,  temperature  histories  at  two  inter- 
nal locations  may  be  used  to  obtain  temperature  or  flux  histories  at 

any  intermediate  points  y,  xi  < y < X2>  by  numerically  solving  the  direct 
problem  on  xi  < x < X2- 

3 Marching  schemes  and  the  inverse  problem 

Let  u’j  denote  u(jAx.nAt).  where  A.r.  Af  are  constant  space  and  time  in- 
crements. and  let  aj  = ({Axy^'{iL'^))/{^t  k{u^)).  3j  = k'{u’J)/k{u'^).  We 
formulate  18  possible  explicit  schemes,  marching  in  the  positive  x direction, 
for  obtaining  u(l,t)  and  Uj(l.t).  The  steady  state  conditions  u{x.O)  = 
Constant.  ?/j(x.O)  = 0.  are  used  throughout.  For  schemes  involving  future 
temperatures,  the  duration  of  the  sensor  record.  Tmax  — -TAf.  is  assumed 
sufficient  to  permit  reconstruction  of  the  surface  waveform  on  a prescribed, 
physically  relevant  time  interval,  Tq  — niAt,  in  < X . For  such  schemes,  it 
is  sometimes  expedient  to  obtain  solution  values  at  or  near  Tmai^  by  simple 
extrapolation  of  immediately  preceding  values.  This  procedure  typically  in- 
duces violent  spurious  oscillations  in  the  vicinity  of  Tmax-  In  our  numerical 
experiments,  provided  the  scheme  was  otherwise  well-behaved,  such  end- 
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point  artifacts  did  not  extend  beyond  5 millisec  to  the  left  of  Tmax^  and 


they  did  not  affect  reconstruction  of  the  surface  history  on  the  remainder 
of  the  interval.  Generally,  the  alternative  of  using  backward  difference  for- 
mulae near  Tnutxi  as  in  method  S5  below,  did  not  prevent  such  end-point 
instabilities.  The  first  three  schemes  below  are  nonhnear  generalizations  of 
methods  presented  in  [1,  Chapter  6]. 

The  first  method,  due  to  D ’Souza,  is  based  on  implicit  time  differencing 
for  the  direct  problem; 
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(3) 


Let  Uj  = Solving  for  we  may  reformulate  (3)  as  a coupled  first 

order  space  marching  system  involving  only  the  nodes  j,  j + 1. 

SI  D’Souza  1. 


u” , n = 1 zV 


Explicit  time  differencing  may  also  be  used  and  leads  to; 
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S2  D’Souza  2. 


= (2-(t; 

)u“  + - u] 

^'j+l 

= 

n — 1 .....  -V  — 1 . 

, .V 

= 2 Ah 

- “hh 

^J+1 

= 2Ah 

“hi 

= g. 

O 

II 

o 

The  'leapfrog'  scheme  is  one  example  of  an  unconditionally  unstable  scheme 
for  the  direct  problem.  It  is  obtained  by  replacing  the  left  hand  side  in  (3) 
with  u^’(  )( - u^“^)/(2At).  As  a space  marching  coupled  system  for 

the  IHCP  this  becomes: 

S3  Leapfrog. 


= 2uy  + i 

^'j+i 

n = 1,...,  A - 1. 

, .V 
^j+i 

= 

Ai 

“A+l 

“j+i 

= U% 

^j+i  - A- 

(6) 


S1-S3  are  0(Ar)  in  the  x-variable  and  lead  to  zz(lA),  u{l  — Ax,t)-  The 
gradient  at  the  wall  can  be  obtained  by  using  backward  differences  in  x. 
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Another  approach  is  to  begin  with  the  heat  conduction  equation  rewritten 
as  a space  marching  first  order  system,  [3],  [7],  and  then  proceed  to  discretize 
that  system.  With  q = k{u)iij:^  we  have 


Ux  = q/k{u 


qj.  = 


Consider  first  O(A.r)  schemes.  The  simplest  example  is  the  following: 
S4  Backward. 


= 

+ Ax  qj 

^+1 

= 

A Ax  Cu'( 

0 

0 

n 

^J+l 

(8) 


A method  due  to  Hills  and  Hensel,  [7],  uses  centered  time  differencing  in  the 
interior  and  backward  time  differencing  at  n — X . 

S5  Hills  and  Hensel,  (1986). 


= li]  + {c7]/2){u^^^  - + A.X-  q]/k{a]) 


:V 


;+i 


C+i 


q]  -f  Ax  - u^-^)/{2Xt), 

uf  + cr^i  uf  - Xx  qf  /k{  Uj 

qf  + A.r 

u^j  + Xx  q^kiu^j),  q^^^  = 

12 


n = 1,.  . .,  iV  - 1, 


(9) 


^'\'e  now  construct  some  further  schemes  based  on  ( < )• 


S6 

Central. 

= u]  + Ax  q^/{k{u]) 

= + )(«"+' 

)/(2Ai).  n = l..... 

.*v  - 1, 

“iVi 

= - '‘jX 

V 

■Ij+i 

= 

- Uy  = Qy 

(10) 

S7 

Future  0. 

^7+1  ~ 

Ax  q’J/ikiu]) 

II 

qj  + Az  + ' - 

- 

- 3ti’])/{-2At).  n = 1, 

. ...  .V  - 2. 

, X-k 

J+l 

2Vi  - Vi 

= 

= 1.0 

^^y  q%l  = Qy 

(11) 

Each  of  S5,  S6,  S7  are  O(Ah^)  schemes  in  the  interior,  with  S7  using  for- 
ward time  differencing.  We  next  consider  0{At)  forward  time  differencing; 
the  forward  formulae  in  S9-S11  are  seldom  used  in  well-posed  problems. 

S8  Future  1. 

Wj+i  = + Sx  q]/k{u^) 
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= 

<,"  + ^x  u.’{u’]){u’]+'  - tt")/At. 

72  = 1.  . . . 

,A’-  1, 

^}+i 

= 

2 “ill’  - 

2<?j+l‘  - -j;+f 

= 

'Jj+i  = ■Jj- 

(12) 

S9 

Future  2. 

ii"  + A,r  q"/k(u") 

^Ij  + i 

5"  + A,c^'(u;)(«"+'^-u")/(2Ai), 

72  = 1,  . 

. .,A-  2, 

r= 

2«;^^-‘  - 

= 

k=l.O 

= 

9?+i  = 9°- 

(13) 

SIO 

Future  3. 

^^J+1 

= 

Uj  + Ax  qjlikiu"-) 

57+ 1 

= 

q]  + Ax  c.'(tt")(u;+^  - u;)/(3At). 

72  = 1,  . 

. .,iV  - 3, 

= 

=: 

29iV-' - fc  = 2,h0 

«“+i  = 9°' 

(14) 
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Sll  Future  4. 

Wj+I  = + Sx  q] /{kiu^j) 


'If+l 

— Q’j  'k  ^’( 

«y(«y-u;)/(-iA(), 

n = 1. . . 

,.,A-  4, 

^J+l 

- _ 
~ -^^+1 

“;  + l 

- 

'/yy-'.  fc  = 3.2.i.o 

«“+i 

II 

►CJ 

o 

+ 

(15) 

To  construct  0(  Aj:^)  schemes,  consider  the  second  order  Taylor  expansion 

= t/"  + Ax(  + 0.5(  Ax)^(f^xx)j 
Q]+i  = + Ax(r^^)y  + 0.5f  A2’)"(c7xx)j  • (16) 

We  may  use  (7)  to  express  a^-x  in  terms  of  Uf,  qt,  as  follows. 

Ihx  = qj:/k{u)  - qk'{ll)u^/P{u) 

= uj(  iL)ut/k{  u)  - q^k'i  u)/k^{  u) 

5xx  — Op'  ( n ^Ux T ^ 

= ^'{u)qut/k{u)  ^ ijj{u)qt/k{u)  - ^'{u)k'{ii)qut/k^{u).  (17) 

Substituting  from  (7)  and  (17),  we  may  replace  the  space  derivatives  on  the 
right  hand  side  of  (16)  with  time  derivatives.  Thus,  with  = Ax/fc(ny), 
b’J  = 0.o{Ax)'^^{u])/k{u]),  c]  = -0.r3{Ax)-k'{u])/k^{u^),  d]  = Axu;{u]), 
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and  e’J  = 0.o{Ax)'^  f^uj'{u'J)/kiu^)  - ^'{u’J)k'{u'J)/k^{u'^)^,  we  obtain 

= it^  + a']q]  + b^{ut)]  + c]{q]f 
Qj^i  = qj+d^iut)]  + b]{qt)]  + e]q]{ut)l  (18) 

Several  schemes  can  be  constructed  from  (18)  depending  on  the  manner  in 
which  the  time  derivatives  (Uf)"  and  {qt)"]  are  approximated.  We  define  7 
such  schemes.  R4  and  R6  - Rll,  by  stipulating  that  the  time  derivatives  in 
( 18)  be  approximated  as  in  S4  and  S6  - Sll  respectively,  with  corresponding 
end-point  conditions. 

4 Lax-Richtmyer  analysis 

The  Lax-Richtmyer  theory  of  difference  . pproximations  is  based  on  Fourier 
analysis  of  the  finearized  problem  with  constant  coefficients,  posed  on  the 
whole  real  t-hne.  In  particular,  the  choice  of  ‘boundary*  approximations  at 
t = 0 and  t = Tmax,  will  play  no  role  whatever  when  this  theory  is  applied  to 
IHCP  space  marching  schemes.  Despite  its  shortcomings,  Fourier  analysis  of 
the  constant  coefficient  model  problem  is  quite  often  a surprisingly  effective 
diagnostic  tool.  Referring  to  equations  (3)  through  (18)  in  Section  3,  let 
3'l  = 0,  aj  = a = (Ar)'V(Af  anz*n),  i-i  - 7 = c"  = 

= 0,  aj  = /i,  = O.oaAt,  and  dy  = yAf.  Following  [12.  Chapters  3.  4], 
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we  now  view  each  space  marching  scheme  as  a 2 x 2 system 


= C{Ax.  j = 0 M - 1, 

^^’0(0  = + -oc  < / < oc.  (19) 


Here.  M is  the  number  of  space  steps  Ax  from  the  interior  sensor  to  the 
active  surface:  C(Ax.  At)  is  a 2 x 2 matrix  of  finite-difference  operators  in 
the  t variable,  assumed  appUed  at  every  point  : /(f)  denotes  exact  initial 
data  for  the  IHCP  such  that  a unique  solution  u{x.  t)  exists  with  sufficiently 
many  bounded  derivatives;  and  S{t),  assumed  small,  represents  the  devia- 
tion of  the  actual  input  data  from  this  exact  data.  Let  Tj{t)  denote  the 

'truncation  error'  on  the  line  /Ax.  and  let  € j{t)  = uq(f)  - u(/Ax.f).  be 

the  difference  between  the  exact  solution  of  the  analytic  problem,  and  the 
numerical  solution  with  noisy  data.  We  have 

6j>i(f)  = C(Ax.  Af)6j(f) -f  Ax  rj(f). 

eo{t)  = Sit).  (20) 

Therefore. 

(j+iit)  = ^(f)  + Ai  y] a-  T,(t). 


(21) 


In  (22),  II  II  denotes  the  norm  in  the  time  variable,  is  the  ‘discrete 
solution  operator’,  and  ||  Tj  ||=  0[:\t  + Ax),  or  smaller,  as  At,  Ax  J.  0. 
Therefore,  for  small  ||  ^ ||,  one  may  expect  reasonable  accuracy  in  the  recon- 
structed surface  histories,  provided  At,  Ax  can  be  chosen  sufficiently  small 
without  making  ||  ||  too  large.  The  latter  cpiantity  measures  the  amount 

by  which  noise  in  the  interior  data  becomes  amplified  when  the  marching 
calculation  reaches  the  active  surface.  We  are  interested  in  comparing  the 
norms  of  discrete  solution  operators  for  the  various  space  marching  schemes 
on  the  same  mesh. 

Fourier  transforming  the  time  variable,  and  putting  9 = <^At,  where  ^ 
is  the  transform  variable,  we  may  find  the  amplification  matrix  G{i\x,9) 
for  each  of  Sl-Sll,  R4,  and  R6-R11.  G(Ax.^)  is  the  Fourier  image  of 
C(Ax,At),  [12,  p.  67],  expressed  in  terms  of  the  normalized  frequency  9. 
Moreover. 

11  11=  max  |G''^(«)|p  = .-I™,,,,  (23) 

0<6I<t 

where  | |/2  denotes  the  Euclidean  norm  of  G.  All  such  matrices  G = [Qij] 
are  2 x 2,  and  we  write  them  in  the  compact  notation  [gw-,  g\2\  921^  922]-  We 
have: 

SI  G = [2  + (T-(7e-'^-l;l,0] 
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52  G = [2-  c^  + (7e'^-l;l,0] 

53  G = [2 lasinO. 

54  G = [l,/z;7(l-e-*^),ll 

55  G = [1  + icfsinO^  /i;  27  sin  9, 1] 

56  G — [1, /i;  Z7.sm^,  1] 

57  G = [Laz;7/;(^)/2,1]  ; p{9)  - - 3 

58  G = [l./i;7(e'^  - 1),1] 

59  G = [l,//;7(e2'^-  l)/2,ll 

510  G = [l.^i;7(e3*^  - l)/3, 1] 

511  G-  [L/^;7(e‘''^-  l)/4,ll 

R4  G = [1  4-  (7{l-  e-‘^)/2,/i;7(l  - + a(l  - e-*^)/2] 

R6  G = [1  + zcr(5in6^)/2, /i;  Z7.szn^,  1 + ^o■(5^7^^^)/2] 

R7  G = [1  + crj[)(^)/4,|r,  7p(^)/2,  1 4-  apiO)!^] 

R8  G = [1  + - l)/2,//;  7(e*^  - 1),  1 + cr(e'^  - l)/2] 

R9  G = [1  + cr(e2»'^  - l)/4,/i;7(e2t^  _ i)/2, 1 -f-  - l)/4] 

RIO  G = [1  4-  cr(e^'^  — l)/6,  7(e^*^  - l)/3, 1 + cr[e^^^  — l)/6] 

Rll  G = [1  4-  - l)/8,//;7(e-^'^  - l)/4, 1 + - l)/8] 

For  the  gun  barrel  IHCP  discussed  in  Sections  1 and  2,  we  have,  in  pre- 
viously indicated  units,  a„un  - 0.004,  At  = 0.1,  and  the  sensor  is  located 
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at  a distance  / = 0.5  from  the  active  surface  2-  = 1.  We  place  M = 1000 
equispaced  nodes  on  /,  so  that  Ax  = 5 x 10“'‘.  In  particular,  this  gives 
Q„„nAt/(A2)^  = 1600,  so  that  the  Courant  condition  for  explicit  time 
marching  schemes  in  the  direct  problem,  [12.  p.  189],  is  severely  violated.  To 
complete  the  parameter  specihcation,  we  give  the  constant  value  3.2  and 
define  knun  = ^'Ctnun  = 0.0128.  With  these  parameter  values,  it  is  instruc- 
tive to  evaluate  at  discrete  points  Oj  = [j  — l)7r/81,  j = 1,82, 

for  each  of  the  above  18  schemes.  This  may  be  done  easily  using  the  soft- 
ware package  MAT  LAB,  [10].  The  results  of  these  computations  are  best 
described  when  broken  up  into  three  groups  as  follows: 

Group  I (Fig.  2A)  O(A^)  Past  Temperature  Schemes:  SI,  S4,  R4. 

Group  II  (Fig.  2B)  0{At^)  Schemes:  S3,  S5,  S6,  S7,  R6,  R7. 

Group  III  (Fig.  2C)  0{i\t)  Future  Temperature  Schemes:  S2,  S8, 

S9,  SIO,  Sll,  R8,  R9,  RIO,  Rll. 

The  Figures  display  To/710  { |G^^°(^)|/2 } versus  the  normalized  frequency 
on  0 < < TT. 

As  is  evident  from  Figs.  2A,  2B,  and  2C,  the  maximum  amplification  fac- 
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tor.  -4njax-  ranges  from  about  10^®  in  Si.  to  about  10'^  in  Sll,  Rll.  More- 
over. only  slight  changes  occur  in  these  Figures  when  the  number  of  space 
nodes  is  reduced  from  1000  to  300.  In  particular,  characteristic  frequency 
domain  signatures  remain  the  same;  Group  I schemes  increase  monotoni- 
cally  with  6,  Group  II  schemes  are  symmetric  about  0 — t/2.  with  a single 
maximum  at  ;r/2.  while  Group  III  schemes  may  have  multiple  maxima.  The 
three  schemes  in  Group  I employ  distinct  space  differencing  methods,  but 
share  backward  time  differencing.  In  Group  II.  S5,  S6,  R6,  have  identi- 
cal traces:  along  with  S3,  these  schemes  use  centered  time  differencing,  but 
each  of  the  four  schemes  uses  a different  space  differencing  technique.  In 
Group  III.  S2  results  from  using  exphcit  differencing  in  D 'Souza's  method 
Si.  This  simple  switch  to  forward  time  differencing  reduces  Amax  by  almost 
ten  orders  of  magnitude,  and  renders  S2  comparable  to  S5,  S6,  R6.  Still 
greater  reductions  are  provided  in  S8,  R8.  by  applying  the  same  switch 
to  S4,  R4.  The  question  arises  as  to  whether  substantial  improvements 
in  S5,  S6,  R6,  might  not  result  from  replacing  centered  time  differencing 
with  a forward  differencing  formula  that  maintains  0(  At^)  accuracy  in  the 
interior.  Schemes  S7,  R7,  were  formulated  with  that  purpose  in  mind. 
However,  as  may  be  seen  from  Fig.  2B.  this  particular  future  temperature 
method  actually  increases  An^ax-  It  is  also  clear  from  the  foregoing  that  the 
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time  differencing  option  plays  a dominant  role  in  determining  frequency  do- 
main signatures,  and  that  some  0{At)  methods  may  be  preferable  to  some 
O(At^)  methods.  The  behavior  of  Si,  S2,  S3.  indicates  that  the  stability 
or  instability  of  these  schemes  in  the  time  marching  direct  problem,  is  of  no 
particular  relevance  in  the  space  marching  inverse  problem. 

The  remaining  six  schemes  in  Group  III.  S9,  R9,  SlO,  RIO,  Sll,  Rll, 
are  based  on  unconventional  forms,  of  forward  differencing  that  induce  cu- 
rious frequency  domain  signatures,  as  well  as  substantial  further  reductions 
in  A,nax-  The  behavior  of  Amax  can  be  correlated  with  the  non  dimen- 
sional parameter  Q = l{ar,unAt)~^^^  mentioned  in  the  Introduction.  We 
have  Q.  = 25  in  Figs.  2A,  2B,  and  2C.  In  Fig.  3A,  Group  I schemes  are 
reexamined  with  / = 1.0  and  all  other  parameter  values  unchanged,  so  that 
A = 50.  In  Fig.  3B,  Group  II  schemes  are  reevaluated  with  At  = 0.01 
and  all  other  parameters  as  in  Fig.  2B.  so  that  Q = 79.  Finally,  in  Fig. 
3C,  Group  III  schemes  are  considered  with  Q = 250,  resulting  from  using 
At  = 0.01,  amtn  = 0.0004,  kmtn  = 0.00128,  and  the  remaining  parameters 
as  in  Fig.  2C.  Together  with  further  computations  involving  various  com- 
binations of  parameters  and  a fair  range  of  Q.  values,  these  Figures  indicate 
an  asymptotic  rule  of  thumb,  A^ax  ~ 10^^,  where  A is  a slowly  varying 
function  of  the  parameters  that  may  be  taken  as  a constant.  Thus,  A % 0.6 


09 


in  Group  I scliemes.  In  Group  II.  A 0.3  for  S5,  S6,  R6.  while  A % 0.35 


for  S7,  R7.  In  Group  HI.  A ^ 0.12  for  Sll,  Rll.  gradually  increasing  to 
about  twice  that  value  for  S8,  R8.  In  Groups  II  and  III.  the  R schemes, 
which  are  second  order  accurate  in  the  space  variable,  are  better  behaved 
than  their  S counterparts.  In  summary,  the  best  schemes  in  Group  II  blow- 
up hke  the  square  root  of  SI.  while  the  best  schemes  in  Group  III  blow-up 
like  the  fifth  root  of  Si. 

5 A numerical  experiment 

It  remains  to  demonstrate  the  relevance  of  the  preceding  analysis  in  the 
computation  of  the  nonlinear  IHCP.  on  a finite  t-interval.  with  end-point 
conditions.  With  reference  to  Fig.  1 in  Section  2.  we  seek  to  recover  rapid 
temperature  and  flux  pulses  at  a:  = 1,  from  highly  attenuated  interior  data. 
Consider  first  the  easier  problem  of  reconstruction  from  = 0.7  when  no 
random  noise  is  added  to  the  data.  In  this  case  Q = 15.  The  sfight  amount 
of  residual  data  noise  originating  from  numerical  computation  of  the  direct 
problem,  is  nevertheless  sufficient  to  trigger  instabilities  if  Group  I schemes 
are  used.  On  the  other  hand.  Group  II  and  III  schemes  perform  weU  un- 
der these  conditions.  For  example,  as  shown  in  Fig.  4,  S5,  S6,  S7,  S8, 
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reproduce  the  exact  temperature  and  gradient  histories.  An  irksome  side 
effect  of  Group  III  schemes  with  multiple  maxima  is  an  overdamping  of  cer- 
tain frequency  ranges,  which  causes  S9,  SlO,  Sll,  to  underestimate  the 
true  solutions.  This  effect  is  more  pronounced  in  the  case  of  the  narrow 
flux  pulse,  than  in  the  case  of  the  temperature  pulse.  For  the  temperature 
history,  the  true  maximum  is  underestimated  by  0.7%  with  S9.  1.1%  with 
SlO,  and  1.4%  with  Sll.  The  corresponding  peak  gradient  underestimates 
are  3.6%).  5.3%,  and  6.8%.  respectively.  The  worst  case,  Sll,  is  shown  in 
Fig.  5.  Note  that  the  last  five  milliseconds  of  each  trace  in  Figs.  4 and  5 
have  been  deleted;  this  is  the  region  where  the  reconstruction  is  seriously 
affected  by  end-point  instabilities. 

The  situation  changes  drastically  when  we  attempt  reconstruction  from 
X = 0.5  with  0.1%  random  noise  added  to  the  interior  data.  Now.  with 
Q.  = 25,  all  Group  I and  Group  II  schemes  are  hopelessly  unstable.  The 
behavior  of  S6,  shown  in  Fig.  6,  is  typical  of  Group  II  schemes  and  leads 
to  amplitudes  on  the  order  of  10®.  In  addition,  most  Group  III  schemes  fail 
to  produce  recognizable  traces  in  this  representative,  yet  difficult,  inverse 
computation.  However,  useful  information  is  recoverable  with  either  of  Sll, 
Rll,  as  may  be  seen  in  Fig.  7.  Although  badly  contaminated  by  noise  and 
unable  to  match  the  true  maximum,  the  reconstructed  gradient  provides 
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a fair  representation  of  the  concentrated  flux  waveform  at  2:  = 1.  The 
slower  temperature  pulse  emerges  relatively  well.  The  result  of  post-filtering 
the  preceding  computation  in  the  time  domain,  using  a three  term  moving 
average,  is  shown  in  Fig.  8.  If  the  level  of  noise  in  the  interior  data  is 
increased  to  1%.  Sll,  Rll  fail  to  produce  satisfactory  results,  even  after 
smoothing.  On  the  other  hand,  with  0.0l9c  added  noise.  S9,  R9  provide 
reconstructions  of  somewhat  better  quality  than  in  Figs.  7 and  S.  and  with 
better  estimates  of  the  true  maxima. 

We  have  not  discussed  regularization  oi  the  IHCP  in  the  present  paper. 
For  the  linear  problem  with  constant  coefficients.  [3].  or.  more  generally, 
whenever  the  analytic  solution  operator  for  the  direct  problem  is  known. 
Tikhonov  regularization  techniques  may  be  applied  to  control  the  growth  of 
errors.  In  the  case  analyzed  in  [3].  such  regularization  is  shown  to  be  equiv- 
alent to  subjecting  the  initial  data  to  a specific  low-pass  filter  in  the  Fourier 
transform  domain.  One  may  also  apply  the  appropriate  fractional  power  of 
that  filter  at  every  step  of  a space  marching  calculation.  For  nonlinear  prob- 
lems. stepwise  filtering  based  on  the  related  linearized  problem  at  each  step 
can  sometimes  be  useful,  [4].  However,  regularization  techniques  in  march- 
ing computations  are  likely  to  be  effective  only  if  the  underlying  difference 
scheme  blows-up  relatively  slowly  with  Q.  The  design  of  such  regularized 
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marching  algorithms  is  another  motivation  for  the  present  study. 
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K / mm 


Figure  1 

Exact  solution  of  the  direct  problem  at  x = 0,  0.5,  1.0. 
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Figure  2 

Behavior  of  on  the  mesh  At  = 0.1,  Ax  = 5 x 10“'^. 
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Figure  3 

Behavior  of  at  various  values  of  H. 
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Figure  4 

Nonlinear  reconstruction  under  favorable  conditions. 
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Figure  5 

Overdaiiiping  in  Sll  leads  to  1.4%  (6.8%)  underestimate  of  peak 
temperature  (gradient)  in  reconstruction  from  x = 0.7. 
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Figure  6 

Behavior  of  S6  under  adverse  conditions  is  typical  of  Group  II  schemes. 
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Figure  7 

Unsmoothed  Sll  reconstruction. 
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Figure  8 

Smootlied  Sll  reconstruction  using  three  term  moving  average. 
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